function [Q,h0Tmp,hxTmp,sigmaTmp] = induceStatVAR(scaling,input)

nx            = size(input.hx,1);
hxTmp         = scaling*input.hx;
eigValues     = eig(hxTmp);
realEigValues = real(eigValues);
imagEigValues = imag(eigValues);
modulus       = sqrt(realEigValues.^2 + imagEigValues.^2);
if any(modulus>=1) || scaling < 0 || scaling > 1
    Q        = 1e6;
    h0Tmp    = NaN(nx,1);
    hxTmp    = NaN(nx,nx);
    sigmaTmp = NaN(nx,nx);
else
   % We have a stationary VAR model
   h0Tmp  = (eye(nx)-hxTmp)*mean(input.factors,2);
   T      = size(input.factors,2);
   residualsTmp = zeros(nx,T);
   for t=2:T
       residualsTmp(:,t) =  input.factors(:,t) - (h0Tmp + hxTmp*input.factors(:,t-1));
   end
   varU     = sum(input.Var_u(:,:,1:T-1),3);
   cov10U   = sum(input.Cov10_u(:,:,2:T),3);
   cov01U   = sum(input.Cov01_u(:,:,2:T),3);
   omegaTmp = 1/(T-1-nx-1)*(residualsTmp*residualsTmp') + 1/(T-1)*(-varU - hxTmp*varU*hxTmp' + cov10U*hxTmp' + hxTmp*cov01U);
   [sigmaTmp,check] = chol(omegaTmp,'lower');
   deScaling     = 1;
   while check ~= 0 
       deScaling = deScaling*0.9;
       omegaTmp = 1/(T-1-nx-1)*(residualsTmp*residualsTmp') + 1/(T-1)*deScaling*(-varU - hxTmp*varU*hxTmp' + cov10U*hxTmp' + hxTmp*cov01U);
       [sigmaTmp,check] = chol(omegaTmp,'lower');
   end
   varX     = reshape((eye(nx^2)-kron3(hxTmp,hxTmp))\reshape(sigmaTmp*sigmaTmp',nx^2,1),nx,nx);
   % Matching all the variances
   Q = 0;
   for i=1:nx
      Q = Q + ((varX(i,i)-input.varXData(i,i))/input.varXData(i,i))^2;
   end
   Q = sqrt(Q/nx)*100;
end


end

